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Obstructive  lung  diseases  in  the  lower  airways  are  a  leading  health  concern  worldwide.  To  improve  our 
understanding  of  the  pathophysiology  of  lower  airways,  we  studied  airflow  characteristics  in  the  lung 
between  the  8th  and  the  14th  generations  using  a  three-dimensional  computational  fluid  dynamics 
model,  where  we  compared  normal  and  obstructed  airways  for  a  range  of  breathing  conditions.  We 
employed  a  novel  technique  based  on  computing  the  Pearson's  correlation  coefficient  to  quantitatively 
characterize  the  differences  in  airflow  patterns  between  the  normal  and  obstructed  airways.  We  found 
that  the  airflow  patterns  demonstrated  clear  differences  between  normal  and  diseased  conditions  for 
high  expiratory  flow  rates  (  >  2300  ml/s),  but  not  for  inspiratory  flow  rates.  Moreover,  airflow  patterns 
subjected  to  filtering  demonstrated  higher  sensitivity  than  airway  resistance  for  differentiating  normal 
and  diseased  conditions.  Further,  we  showed  that  wall  shear  stresses  were  not  only  dependent  on 
breathing  rates,  but  also  on  the  distribution  of  the  obstructed  sites  in  the  lung:  for  the  same  degree  of 
obstruction  and  breathing  rate,  we  observed  as  much  as  two-fold  differences  in  shear  stresses. 
In  contrast  to  previous  studies  that  suggest  increased  wall  shear  stress  due  to  obstructions  as  a  possible 
damage  mechanism  for  small  airways,  our  model  demonstrated  that  for  flow  rates  corresponding  to 
heavy  activities,  the  wall  shear  stress  in  both  normal  and  obstructed  airways  was  <  0.3  Pa,  which  is 
within  the  physiological  limit  needed  to  promote  respiratory  defense  mechanisms.  In  summary,  our 
model  enables  the  study  of  airflow  characteristics  that  may  be  impractical  to  assess  experimentally. 
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1.  Introduction 

Obstructive  lung  diseases  are  one  of  the  most  common  and  life- 
threatening  conditions,  causing  high  morbidity  and  mortality  in 
industrialized  countries  [1],  Many  lung  diseases,  such  as  asthma, 
bronchitis,  and  chronic  obstructive  pulmonary  disease  (COPD),  are 
characterized  by  the  obstruction  of  airflow  in  the  lungs.  According 
to  the  World  Health  Organization,  COPD  is  the  fourth  leading 
cause  of  death  worldwide  and  is  estimated  to  become  the  third 
leading  cause  by  2030  [f  ],  The  dominant  cause  of  such  obstructive 
lung  diseases  is  exposure  to  lung  irritants,  resulting  in  the 
contraction  of  smooth  muscle  in  airway  walls  and  inflammation, 
which  in  turn  lead  to  narrowing  of  the  airway  lumen  [2], 
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The  development  of  chronic  disease  conditions  associated  with 
airway  obstruction  may  be  altered  by  proper  and  timely  inteiwen- 
tion  [3].  Thus,  early  detection  and  correct  diagnosis  are  critical. 
However,  some  airway  conditions  remain  hard  to  diagnose, 
especially  those  pertaining  to  the  lower  airways  [4],  For  example, 
significant  discordance  has  been  reported  between  commonly 
used  references  for  the  interpretation  of  spirometry  data  [5],  and 
the  American  Thoracic  Society  does  not  recommend  the  use  of 
spirometric  parameters,  such  as  the  mid-expiratory  flow,  to  detect 
lower  airway  disease  [6],  According  to  the  Third  National  Health 
and  Nutrition  Examination  Survey  (NHANES  111),  undiagnosed 
airflow  obstruction  (12%)  is  more  common  than  diagnosed  COPD 
(3.1%)  or  asthma  (2.7%)  in  the  United  States  (US.)  and  is  associated 
with  impaired  health  and  functional  status  [4].  Undiagnosed 
pulmonary  conditions  have  also  been  reported  in  U.S.  military 
populations  [7].  Out  of  105  active  duty  military  patients  who 
complained  of  exertional  dyspnea,  one  quarter  (25  patients)  had 
no  specific  diagnosis  even  after  a  comprehensive  standard  evalua¬ 
tion  [7], 

Understanding  the  airflow  characteristics  in  the  lower  air¬ 
ways  of  the  lungs  could  provide  detailed  insights  on  respiratory 
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physiology  and  pathophysiology  and  be  potentially  useful  for 
diagnosis  of  lower  airway  diseases.  However,  due  to  limitations 
in  spatial  and  temporal  resolutions,  it  is  difficult  to  study  local 
airflow  characteristics  in  the  lower  airways  using  experimental 
techniques.  For  example,  a  pulmonary  function  test  (PFT)  involves 
measuring  the  volume  of  airflow  during  forced  expiration  at  the 
mouth  and  comparing  it  with  standardized  normal  values  for 
diagnosis.  However,  because  PFT  measures  airflow  in  the  bulk 
phase,  it  is  insensitive  to  regional  differences  in  airflow  patterns, 
especially  in  the  lower  airways.  Similarly,  although  hyperpolarized 
noble  gas  magnetic  resonance  imaging  can  provide  static  and 
dynamic  ventilation  maps  of  the  lungs  [8],  it  is  still  not  possible  to 
obtain  details  of  airflow  characteristics,  such  as  velocity  contours, 
in  the  lower  airways  [9]. 

Computational  fluid  dynamics  (CFD)  provides  the  means  to 
investigate,  in  detail  and  under  controlled  conditions,  airflow 
characteristics  in  three-dimensional  (3-D)  models  of  the  respira¬ 
tory  system.  Previous  CFD  modeling  efforts  in  pulmonary  health 
have  focused  mainly  on  (1)  establishing  particle  dosimetry  in 
normal  airways  to  optimize  inhaled  drug  delivery  and  understand 
risks  due  to  airborne  pollutants  [10,11],  (2)  understanding  sleep 
apnea  associated  with  obstructions  in  extrathoracic  airways 
[12,13],  and  (3)  investigating  airway  obstructions  in  the  upper 
generations  of  the  lungs  [14,15].  In  contrast,  we  focus  on  exploring 
airflow  characteristics  in  the  lower  airways  under  different  breath¬ 
ing  conditions  and  examining  influences  of  lower-airway  obstruc¬ 
tions  as  compared  with  normal  conditions.  Moreover,  we  aim  to 
characterize  airflow  during  both  the  inhalation  and  exhalation 
phases  of  respiration,  as  opposed  to  previous  studies  which  have 
mostly  investigated  the  inhalation  phase  [16,17].  Specifically,  we 
developed  CFD  models  of  normal  and  obstructed  airways  between 
the  8th  and  the  14th  generations  with  different  obstruction 
patterns.  We  then  compared  the  airflow  patterns  and  flow- 
induced  wall  shear  stresses  for  different  flow  conditions.  Our 
results  quantify  differences  in  airflow  characteristics  due  to 
obstructions  commonly  associated  with  lower  airway  lung  dis¬ 
eases  for  both  phases  of  respiration. 


2.  Methods 

We  developed  two  sets  of  airway  models  to  validate  our 
approach  and  study  airflow  characteristics.  The  first  set  of  models 
(central  and  peripheral  airways),  which  included  airway  genera¬ 
tions  from  the  trachea  down  to  the  14th  generation,  was  used  to 
validate  the  computational  approach  of  using  a  reduced  number  of 
airway  generations  to  model  airflow  in  the  lungs.  The  second  set  of 
models  (normal  and  obstructed  airways)  included  only  the  airway 
branches  between  the  8th  and  14th  generations  and  was  used  to 
study  the  effects  of  obstructions  in  these  lower  airways  on  the 
airflow  characteristics.  To  quantify  the  similarity  and  dissimilarity 
in  airflow  patterns  in  normal  and  obstructed  airways,  we  used 
correlation  coefficients  and  filtering  techniques.  We  used  different 
geometries  and  boundary  conditions  for  these  two  sets  of  models, 
which  are  described  in  detail  below. 

2.1.  Airflow  models  of  central  and  peripheral  airways 

To  assess  the  validity  of  using  airway  models  with  a  reduced 
number  of  airway  branches,  we  developed  three  3-D  airway 
models  for  central  and  peripheral  airways  with  different  numbers 
of  airway  branches.  Fig.  1  shows  the  bronchi  model  (A),  the  central 
airway  model  (B),  the  central-peripheral  airway  model  (C),  and 
airway  geometiy  and  experimental  data  for  in  vitro  velocity 
measurements  by  Isabey  and  Chang  [18]  (D).  The  bronchi  model 
mimics  the  physical  model  from  Isabey  and  Chang  [18],  which  was 


B 


Central  airway  model 


Central-peripheral  airway  model 
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Fig.  1.  Computational  models  of  the  bronchi  (A),  central  (B),  and  central-peripheral 
airways  (C)  as  well  as  experimental  data  (D)  from  Isabey  and  Chang  [18].  Airway 
geometries  for  both  computational  and  experimental  models  were  based  on  the 
Horsfield  model  [19,20].  Left:  pressure  fields  obtained  from  computational  fluid 
dynamics  (CFD)  simulations  (A)-(C)  and  airway  geometry  used  in  the  experiment 
(D).  Right:  contours  of  axial  velocity  normalized  by  the  maximum  at  the  cross- 
section  (a-a')  obtained  from  CFD  simulations  (A)-(C)  and  the  experiment  (D). 
In  the  plots,  each  consecutive  contour  represents  a  velocity  increment  of  0.2.  (For 
interpretation  of  the  references  to  color  in  this  figure  legend,  the  reader  is  referred 
to  the  web  version  of  this  article.) 


used  to  measure  velocity  profiles,  and  consists  of  branches  from 
the  trachea  (i.e.,  the  1st  generation)  to  five  lobular  bronchi.  Isabey 
and  Chang  [18]  used  a  truncated  representation  of  the  Horsfield 
model  [19,20],  which  accurately  describes  healthy  lung  anatomy 
obtained  from  a  lung  cast  of  a  male  cadaver  and  captures  the 
asymmetry  of  the  lung  structure.  For  the  central  airway  model 
(Fig.  1(B)),  we  extended  the  bronchi  model  by  including  airway 
branches  from  the  trachea  to  the  7th  airway  generation  to 
replicate  the  complete  central  ainway  model  available  from  Hors- 
fleld  [19,20]  (see  Tables  SI  and  S2  and  Fig.  SI  in  Supplementary 
material).  The  central-peripheral  airway  model  represents  both 
the  central  airways  and  one  set  of  peripheral  airways.  The 
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peripheral  airway  component  starts  at  one  branch  of  the  7th 
generation  in  the  central  airway  model  and  includes  all  airway 
branches  between  the  8th  and  14th  generations  (Fig.  1(C)).  Due  to 
the  unavailability  of  data  regarding  the  branching  angles  in  the 
lower  airways,  the  peripheral  airways  were  assumed  to  be  sym¬ 
metric  [21,22].  Moreover,  the  effect  of  asymmetry  on  flow  patterns 
was  shown  to  be  insignificant  for  low  Reynolds  number  flow 
associated  with  the  peripheral  airways  [23]. 

We  created  the  3-D  airway  geometries  using  computer-aided 
design  software  (AutoCAD  2012,  The  Autodesk).  We  used  airway 
dimensions,  such  as  diameter,  length,  branching  angle,  and 
branching  radius  (curvature),  for  each  branch  in  the  central 
airways  [20]  and  average  values  of  airway  dimensions  for  each 
generation  in  the  peripheral  airways  [19]  (see  Tables  SI  and  S2  and 


Fig.  SI  in  Supplementary  material).  Transitional  zones  near  the 
bifurcations  were  built  by  lofting  circular  cross-sections  of  mother 
and  daughter  branches  in  AutoCAD.  Due  to  the  complexity  of  the 
geometries  obtained  after  lofting,  we  did  not  incorporate  the 
rounding  of  the  carinal  ridges  [24[. 

In  this  set  of  simulations,  our  aim  was  to  compare  the  airflow 
patterns  predicted  by  our  models  with  in  vitro  experimental  data 
for  steady-state  inspiratory  flow  from  Isabey  and  Chang  [18[. 
To  mimic  the  experimental  conditions  [18],  we  applied  a  velocity 
boundary  condition  at  the  inlet  of  the  1st  generation  (trachea)  and 
outflow  boundary  conditions  at  the  end  of  the  terminal  branches 
of  the  airway.  A  uniform  velocity  magnitude  (1  m/s)  normal  to  the 
boundary  was  specified  at  the  inlet  emulating  the  experiments. 
This  corresponded  to  a  flow  rate  of  200  ml/s,  which  is  close  to  the 
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Fig.  2.  Three-dimensional  (3-D)  geometries  of  peripheral  airways  (i.e.,  8th  generation  to  14th  generation).  (A)-(E);  detailed  description  of  normal  and  obstructed  branches  at 
the  nth  airway  generation  (A)  and  3-D  geometries  of  normal  airway  (B),  symmetric  obstruction  (C),  asymmetric  obstruction  (D),  and  random  obstruction  (E)  models. 
dn,  branch  diameter;  Rn,  bifurcation  diameter;  0,  branching  angle.  The  dimensions  of  all  normal  branches  were  taken  from  the  Horsfield  model  [19].  Values  for  0  ( =45  )  and 
Rn  (=2dn)  were  chosen  within  the  physiological  range  [19].  The  diameter  of  the  obstructed  branch  was  reduced  by  50%.  In  the  normal  airway  model,  all  branches  were 
unobstructed.  In  the  symmetric  obstruction  model,  all  branches  were  obstructed.  In  the  asymmetric  obstruction  model,  one  of  the  second  generation  branches  and  all  of  its 
following  daughter  branches  were  obstructed.  In  the  random  obstruction  model,  obstruction  sites  were  chosen  randomly.  In  all  models,  airway  branches  in  the  first  and  last 
generations  were  unobstructed. 
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mean  inspiratoiy  flow  rate  in  sedentary  condition  [25],  We  did  not 
add  entrance  lengths  to  the  airway  models  because  an  extended 
entrance  length  was  not  included  in  the  experiments.  Flow  was 
assumed  to  distribute  unequally  at  the  outlets  to  emulate  the 
experiments  and,  therefore,  flow  rate  weighting  was  specified  as  a 
constraint  in  the  outflow  boundary  condition  at  each  of  the 
terminal  branches.  We  estimated  the  flow  rate  weighting  based 
on  the  volume  fraction  of  the  lung  supplied  by  the  terminal 
branch,  as  done  in  the  experiment  [18],  First,  we  determined  the 
flow  rate  weighting  in  each  of  the  five  lobar  bronchi  using  their 
corresponding  lobar  volume  fractions  from  the  literature  [20]:  20% 
(upper  right),  10%  (middle  right),  25%  (lower  right),  20%  (upper 
left),  and  25%  (lower  left).  Next,  we  divided  the  flow  in  each  lobar 
bronchus  by  1/2"  to  obtain  flow  rates  at  the  outlets  in  the 
nth  generation.  We  prescribed  a  no-slip  boundary  condition  at 
the  walls. 

For  numerical  calculations,  we  created  unstructured  meshes 
using  GAMBIT  2.4.6  (ANSYS  Inc.,  Canonsburg,  PA).  We  used 
2  million,  4  million,  and  13  million  cells  for  the  bronchi,  central 
airway,  and  central-peripheral  aiiway  geometries,  respectively, 
based  on  grid-independence  tests  (see  Fig.  S2  and  Table  S3  in 
Supplementary  material).  To  obtain  a  numerical  solution  for  the 
Navier-Stokes  equations  governing  airflow,  we  used  the  steady- 
state,  laminar  flow  solver  in  FLUENT  14  (ANSYS  Inc.,  Canonsburg, 
PA).  The  assumption  of  laminar  flow  was  valid,  as  the  highest 
Reynolds  number  for  the  flow  was  1055  at  the  1st  generation;  it 
was  lower  for  subsequent  airway  generations.  For  each  simulation, 
we  used  8  2.8-GHz  Intel  Xeon  quad-core  Nehalem  processors  at 
the  Department  of  Defense  (DoD)  Supercomputing  Resource 
Center  located  at  the  U.S.  Army  Research  Laboratory. 


2.2.  Airflow  models  for  normal  and  obstructed  airways 

We  examined  airflow  patterns  in  healthy  and  diseased  airways 
by  developing  corresponding  3-D  models  for  normal  (unob¬ 
structed)  and  obstructed  airways.  Because  airway  obstruction 
occurs  mainly  in  the  peripheral  airways  for  most  COPD 
cases  [26],  we  created  peripheral  airway  geometries  from  the 
8th  to  14th  generation  [20]  without  truncation  of  airway  branches. 
To  examine  the  effects  of  disease  location  and  severity  on  airflow 
patterns,  we  constructed  one  normal  and  three  different 
obstructed  airway  geometries,  consisting  of  symmetric,  asym¬ 
metric,  and  random  obstructions.  Fig.  2  shows  the  geometric 
configurations  representing  normal  and  obstructed  branches  (A), 
and  the  geometries  of  the  normal  airways  (B),  the  symmetric 
obstruction  (C),  the  asymmetric  obstruction  (D),  and  the  random 
obstruction  (E).  In  the  normal  airway  geometry,  none  of  the  airway 
branches  was  obstructed.  In  the  obstructed  airway  geometries,  we 
constricted  some  of  the  branches  between  the  9th  and  13th 
generations,  which  have  been  identified  by  a  bronchoscopy  study 
[27]  as  the  most  frequent  disease  sites.  In  the  symmetric  obstruc¬ 
tion  model,  each  of  the  64  airway  branches  between  the  9th  and 
13th  generations  was  obstructed.  In  the  asymmetric  obstruction 
model,  only  the  9th  generation  airway  branch  on  the  left  (i.e.,  in 
the  +x  direction)  and  its  daughter  branches  were  obstructed, 
resulting  in  a  total  of  32  obstructions.  In  the  random  obstruction 
model,  we  randomly  constricted  32  airway  branches  of  the  64 
branches  between  the  9th  and  13th  generation  (Eig.  2(E)).  Eigh¬ 
teen  branches  on  the  left  (i.e.,  in  the  +x  direction)  and  14  branches 
on  the  right  (i.e.,  -x  direction)  were  obstructed  (see  Fig.  S3  in 
Supplementary  material).  The  resulting  volume  reduction  in  our 
geometries  was  17%  for  the  symmetric  obstruction  model  and  8% 
for  the  asymmetric  and  random  obstruction  models  (Table  1 ).  This 
is  consistent  with  airway  surface  area  reductions  in  disease 
conditions  reported  by  histological  studies  ]26,28]. 


In  this  set  of  simulations,  we  used  pressure  boundaiy  condi¬ 
tions  at  the  inlet  and  outlets  as  opposed  to  the  velocity  and 
outflow  boundary  conditions  used  for  the  first  set  of  simulations 
(i.e.,  for  central  and  peripheral  airways).  By  using  the  pressure 
boundary  conditions,  we  could  compute  the  change  in  air  flow 
rate  due  to  obstructions  in  the  lower  airways  by  subjecting  the 
normal  and  obstructed  airways  to  the  same  pressure  differential. 
An  implicit  assumption  was  that  both  heaithy  and  diseased  lungs 
would  experience  the  same  pressure  differential  across  the  lungs 
for  a  given  breathing  effort.  Within  the  physiological  range  [29], 
we  applied  pressure  differentials  (AP)  of  18, 10,  and  2  Pa  to  each 
airway  geometry  for  both  inspiratory  and  expiratory  flows  to 
simulate  different  breathing  conditions.  Table  1  shows  the  result¬ 
ing  flow  rates  at  the  inlet  of  peripheral  airways  and  the  estimated 
flow  rates  at  the  trachea,  which  correspond  to  a  wide  range  of 
breathing  conditions  between  sedentary  (~  300  ml/s)  to  forced 
breathing  (~5000  ml/s)  ]30,31].  To  estimate  the  flow  rates  at  the 
trachea,  we  assumed  that  similar  sets  of  peripheral  airways  follow 
each  of  the  central  airway  outlets  as  has  been  used  in  previous 
airway  models  [32,33],  Accordingly,  we  estimated  flow  rates  at  the 
trachea  to  be  2^  times  the  values  at  the  peripheral  airways.  We 
extended  the  inlet  of  the  peripheral  airway  to  impose  a  fully 
developed  flow.  To  this  end,  we  added  an  entrance  length  of 
300  mm  (~61  x  diameter)  to  the  first  branch  of  the  peripheral 
airway  (i.e.,  the  8th  generation  from  the  trachea).  The  minimum 
entrance  length  required  to  obtain  a  fully  developed  flow  based  on 
the  maximum  Reynolds  number  was  ~166  mm.  Because  we  used 
pressure  boundary  conditions,  the  flow  velocities  were  not  known 
a  priori  and,  therefore,  we  included  an  extended  entrance  length  to 
ensure  fully  developed  flow  for  a  range  of  applied  pressure 
differentials.  We  also  applied  a  no-slip  boundary  condition  at 
the  walis. 

We  created  unstructured  meshes  with  ~  7  million  cells  for  each 
of  the  airway  geometries,  based  on  grid-independence  tests  (see 
Fig.  S2  and  Table  S3  in  Supplementary  material).  The  flows  in  our 
models  were  laminar  with  a  maximum  Reynolds  number  of  563  at 
the  1  St  generation  of  the  peripheral  airway,  corresponding  to  the 
peak  inspiratory  flow  rate  of  32.7  ml/s  (Table  1)  and  maximum 
Reynolds  numbers  of  374,  201,  and  241  for  the  symmetric,  asym¬ 
metric,  and  random  obstruction  models,  respectively.  The  max¬ 
imum  Reynolds  number  for  the  obstructed  airway  models  was 
observed  at  the  constrictions  at  the  uppermost  generation  (i.e.,  the 
9th  generation  for  the  symmetric  and  asymmetric  obstruction 
models  and  at  the  10th  generation  for  the  random  obstruction 
model).  A  steady-state  solver  was  used  with  a  second-order 
upwind  scheme  for  the  convective  terms  and  the  SIMPLEC  algo¬ 
rithm  for  pressure-velocity  coupling.  Each  simulation  performed 
at  the  DoD's  Supercomputing  Resource  Center  took  approximately 
2-3  days. 

2.3.  Correlation  coefficient  as  a  similarity  measure 

We  quantified  differences  between  airflows  in  the  obstructed 
and  the  normal  airway  models  by  computing  Pearson's  correlation 
coefficient  for  velocity  contours  at  a  particular  cross-sectional 
plane.  For  numerical  calculations,  we  obtained  two  matrices, 
A  and  B,  where  matrix  elements  represent  the  velocity  magnitudes 
at  each  x-y  spatial  coordinate.  We  discretized  the  cross-sectional 
plane  in  the  region  of  interest  by  creating  fine  meshes.  At  each  grid 
point  (Xj,  yj),  we  assigned  the  corresponding  velocity  magnitude 
(v„  Vj)  as  a  matrix  element  (i,  j).  The  correlation  coefficient  (r)  was 
defined  [34]  as  follows: 

liZiiAi-AKBii-B) 

r  =  ,- . -  -  - - j-jg,  ( 1 ) 

^JZi'LiiAij-Af'LiljlBij-Bf 
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Table  1 

Airway  volume,  flow  rate,  and  resistance  of  the  peripheral  airways  (i.e.,  8th-14th  generation). 

Normal  Symmetric  Asymmetric  Random 

_ Obstruction _ Obstruction _ Obstruction 

Airway  voiume,  ml  (ratio  to  the  normal  geometry,  %) 

10.6  (100)  8.8  (83)  9.7  (92)  9.7  (92) 


Flow  rate,  ml/s  (flow  rate  at  the  trachea,*  ml/s) 


Inspiratory  flow, 
AP  (Pa) 

18 

10 

2 

32.7  (4,200) 

19.6  (2,500) 

4.3  (550) 

18.6  (2,400) 

11.8  (1,500) 

3.0  (380) 

29.5  (3,800) 

17.7  (2,300) 

4.0  (510) 

27.3  (3,500) 

16.5  (2,100) 

3.8  (490) 

Expiratory  flow, 
AP  (Pa) 

18 

29.3  (3,800) 

18.0  (2,300) 

26.9  (3,400) 

25.8  (3,300) 

10 

18.2  (2,300) 

11.7  (1,500) 

16.7  (2,100) 

16.0  (2,000) 

2 

4.1  (520) 

3.0  (380) 

4.0  (510) 

3.8  (490) 

Airway  resistance,  craH20-s/l  (ratio  to  the  normai,  %) 

Inspiratory  flow, 
AP  (Pa) 

18 

5.61  (100) 

9.86  (176) 

6.22  (111) 

6.73  (120) 

10 

5.22  (100) 

8.61  (165) 

5.77  (111) 

6.16  (118) 

2 

4.72  (100) 

6.78  (144) 

5.17  (110) 

5.42  (115) 

Expiratory  flow, 
AP  (Pa) 

18 

6.26  (100) 

10.17  (163) 

6.83  (109) 

7.13  (114) 

10 

5.59  (100) 

8.75  (156) 

6.09  (109) 

6.38  (114) 

2 

4.93  (100) 

6.78  (138) 

5.17  (105) 

5.42  (110) 

AP,  change  in  pressure.  *Estimated  by  multiplying  the  flow  rate  of  peripheral  airways  by  2’. 


where  A  and  B  denote  mean  velocity  values  of  matrices  A  and  B, 
respectively.  We  used  ~  30,000  grid  points  for  each  matrix  based 
on  grid-independence  tests  of  correlation  coefficient  values. 


3.  Results 

3.1.  Inspiratory  airflow  in  central  and  peripheral  airways 

Modeling  the  entire  airway  geometry  is  challenging  both 
experimentally  and  computationally.  The  human  respiratory  sys¬ 
tem  contains  >10^  conducting  ainvay  branches  with  ~  20-fold 
size  differences  between  the  trachea  and  terminal  bronchioles 
[20,32],  Thus,  the  human  respiratory  system  is  often  studied  using 
airway  geometries  with  a  reduced  number  of  branches  [33,35], 
Even  though  Choi  et  al,  [36]  examined  the  effects  of  upper  airway 
truncation  on  downstream  airflow,  it  is  not  well  established  how 
modeling  a  reduced  number  of  outlet  branches  affects  the  down¬ 
stream  airflows.  To  address  this  issue,  we  simulated  airflows  using 
the  bronchi,  central,  and  central-peripheral  airway  models  and 
compared  their  flow  patterns. 

We  examined  the  flow  in  the  left  main  bronchus,  which  is  not 
associated  with  the  outlets  in  any  of  the  models,  to  avoid 
immediate  boundary  effects.  Fig,  1  shows  pressure  and  velocity 
fields  of  the  bronchi,  central  airway,  and  central-peripheral  airway 


models.  For  comparisons  with  the  experimental  data  from  Isabey 
and  Chang  [18],  we  followed  their  method  of  normalizing  the 
contour  plots  of  axial  velocity  with  respect  to  mean  values.  In  all 
models,  pressure  drops  across  the  same  airway  branches  were  not 
significantly  different.  Similarly,  the  velocity  fields  at  the  main 
bronchus  did  not  differ  considerably  across  models.  In  all  models, 
the  velocity  magnitude  approached  zero  near  the  wall,  as  expected 
from  the  no-slip  boundaiy  condition.  Also,  it  was  commonly 
shown  in  all  three  models  that  the  high-velocity  domain  was 
skewed  toward  the  lower  left  bronchus,  which  was  similarly 
observed  in  a  previous  computational  airflow  study  [24]  for 
non-planar  Horsfield  model.  All  models  especially  displayed 
horseshoe-shaped  velocity  contours  from  1,2  to  1,6  at  the  left 
bronchus,  with  minor  differences  in  their  width  across  the  a-a' 
axis.  This  suggests  that  reduced  branch  geometries  with  appro¬ 
priate  boundary  conditions  are  acceptable  in  modeling  lung  air¬ 
ways  if  the  region  of  interest  is  not  at  the  outlet  branches. 

We  validated  our  models  by  comparing  the  velocity  contours  at 
the  left  main  bronchus  with  experimental  data.  Fig,  1(D)  shows  the 
geometry  used  in  an  experiment  by  Isabey  and  Chang  [18]  to 
measure  velocity  profiles  in  an  in  vitro  model  and  the  experi¬ 
mental  data  of  normalized  axial  velocity  contours.  The  contours 
from  the  models  and  the  experiment  had  similar  patterns  with 
minor  discrepancies.  As  shown  in  Fig,  1(D),  the  high-velocity 
domain  was  slightly  narrower  across  the  a-a'  axis  compared  with 
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Fig.  3.  High  inspiratory  steady  flow  in  the  normal  airway  (A),  symmetric  obstruc¬ 
tion  (B),  asymmetric  obstruction  (C),  and  random  obstruction  (D)  models.  The 
applied  pressure  drop  across  the  extended  geometry  was  18  Pa.  Pressure  fields  and 
velocity  fields  of  the  coronal  and  uppermost  axial  planes  are  shown  on  the  left  and 
right,  respectively.  (For  interpretation  of  the  references  to  color  in  this  figure 
legend,  the  reader  is  referred  to  the  web  version  of  this  article.) 


Fig.  4.  High  expiratory  steady  flow  in  the  normal  airway  (A),  symmetric  obstruction  (B), 
asymmetric  (C),  and  random  obstruction  (D)  models.  Tire  applied  pressure  drop  across 
the  extended  geometry  was  18  Pa.  Pressure  fields  and  velocity  fields  of  the  coronal  and 
uppermost  axial  planes  are  shown  on  the  left  and  right,  respectively.  Centers  of  vortices 
are  marked  by  “+.”  (For  interpretation  of  the  references  to  color  in  this  figure  legend, 
the  reader  is  referred  to  the  web  version  of  this  article.) 


those  predicted  by  our  models  (Fig.  1(A)-(C)).  We  attributed  these 
discrepancies  to  the  attachment  of  the  linear  resistors  in  the 
experimental  setup,  an  artificial  effect  not  explicitly  accounted 
for  in  the  model.  Moreover,  the  airflow  patterns  in  the  experi¬ 
mental  data  did  not  reflect  the  symmetry  in  the  airway  geometry 
about  the  a-a'  axis  [as  assumed  in  the  experimental  model  [18]], 
which  suggested  limited  precision  in  the  experimental  data.  Some 
of  the  differences  between  the  model  and  experimental  airflow 
patterns,  therefore,  may  be  attributed  to  experimental 


discrepancies.  However,  overall,  we  found  that  our  models  and 
the  experiment  exhibited  a  good  qualitative  agreement  between 
flow  patterns. 

3.2.  Airflow  in  normal  and  obstructed  airways 

Airway  resistance  is  a  measure  of  the  opposition  to  the  airflow 
caused  by  geometric  properties,  such  as  airway  obstruction,  and  is 
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used  to  characterize  and  evaluate  obstmctive  lung  diseases.  It  is 
defined  as  the  ratio  of  the  total  pressure  drop  to  the  flow  rate  [37], 
To  calculate  peripheral  airway  resistance,  we  included  only  the 
peripheral  ain/vay  geometry  (i.e.,  the  region  below  b-b'  in 
Figs.  3  and  4)  and  extended  the  inlet  of  the  first  branch  to  obtain  a 
fully  developed  flow.  Table  1  shows  the  flow  rate  and  resistance  values 
of  each  airway  model  at  different  applied  pressure  drops.  Resistance 
values  were  dependent  on  the  degree  and  geometric  distribution  of 
the  obstruction  sites.  In  the  symmetric  obstruction  model,  resistance 
values  were  highest  among  all  the  models,  as  expected  fi'om  the 
largest  reduction  in  airway  volume.  The  resistance  values  for  the 
random  and  asymmetric  obstruction  models  were  close  to  each  other 
as  it  should  be  expected  from  the  similar  amount  of  airway  volume 
reduction  (Table  1 );  nevertheless,  the  values  for  the  random  obstruc¬ 
tion  model  were  slightly  higher  (4-8%)  than  the  asymmetric  obstruc¬ 
tion  values  depending  on  the  pressure  drop  applied.  It  should  be  noted 
here  that  only  the  distribution  of  the  obstructions  was  different  in  the 
asymmetric  and  random  obstruction  models:  in  the  asymmetric 
obstruction  model,  all  32  branches  on  the  left  were  obstmcted, 
whereas  for  the  random  obstruction  model,  18  branches  on  the  left 
and  14  branches  on  the  right  were  obstructed. 

Resistance  values  increased  with  higher  pressure  drop  for  both 
inspiratory  and  expiratory  flows,  consistent  with  a  previous 
model  [38],  In  addition,  the  difference  in  resistance  values 
between  the  normal  and  obstructed  airway  models  increased  with 
higher  flow  rates.  For  the  high  inspiratory  flow  driven  by  an  18-Pa 
pressure  drop,  resistance  values  for  the  symmetric,  asymmetric, 
and  random  obstruction  models  were  76%,  11%,  and  20%  higher 
than  for  the  normal  airway  model,  respectively  (Table  1 ).  For  the 
low  inspiratory  flow  driven  by  a  2-Pa  pressure  drop,  resistance 
values  of  the  symmetric,  asymmetric,  and  random  obstruction 
models  were  44%,  10%,  and  15%  higher  than  the  normal  airway 
model,  respectively  (Table  1).  Similarly,  for  expiratory  flows, 
differences  between  the  normal  and  obstructed  airway  models 
increased  with  flow  rate  (Table  1 ). 

Increases  in  the  resistance  values  for  higher  pressure  drops, 
corresponding  to  a  higher  flow  rate,  were  due  to  increases  in  the 
secondary  flow.  Based  on  the  Hagen-Poiseuille  equation,  which 
assumes  a  fully-developed  axial  flow  condition,  the  airway  resis¬ 
tance  should  be  constant  for  the  same  geometry.  However, 
because  the  airway  geometries  in  our  models  constitute  a  3-D 
network  of  conduits  that  are  curved  near  the  bifurcations,  sec¬ 
ondary  flows  were  expected  to  develop,  as  reported  in  the 
literature  [39],  In  such  a  system,  flow  resistance  does  not  remain 
constant  (as  expected  from  the  Hagen-Poiseuille  equation)  as  the 
secondary  flow  increases  with  an  increase  in  the  pressure  drop. 
Moreover,  the  effect  of  the  secondary  flow  is  more  prominent  near 
the  constrictions  in  the  obstructed  airways  due  to  the  abrupt 
changes  in  the  geometry.  This  explains  the  increase  of  differences 
in  the  airway  resistance  between  the  obstructed  and  normal 
airways  with  an  increase  in  the  pressure  drop. 

Fig.  3  shows  pressure  and  velocity  fields  for  high  inspiratory 
flow  driven  by  a  total  pressure  drop  of  18  Pa.  Here,  the  airway 
models  consist  of  the  peripheral  airways  as  well  as  the  extended 
entrance  lengths  used  to  obtain  a  fully  developed  flow  profile. 
As  resistance  in  the  peripheral  airways  increased  due  to  obstruc¬ 
tion  compared  with  the  normal  airway  model,  we  observed  higher 
pressure  gradients  in  the  peripheral  airways  (Fig.  3).  Consequently, 
the  pressure  drop  in  the  extended  entrance  sections  (not  shown) 
of  the  airways  decreased  to  compensate  for  the  increase  in 
pressure  drop  across  the  peripheral  airways.  The  resulting  pres¬ 
sure  distribution  across  the  peripheral  airways  was,  therefore, 
different  between  the  models  even  though  the  same  total  pressure 
boundary  condition  of  18  Pa  was  applied  to  all  models. 

We  then  focused  on  airflow  patterns  to  gauge  how  obstructions, 
modeled  as  geometric  changes  in  the  lower-generation  branches. 


manifested  themselves  as  altered  flow  in  an  8th  to  14th  generation 
airway  model.  Specifically,  we  compared  velocity  contours  at  a 
coronal  (i.e.,  perpendicular  to  the  z-axis)  and  an  axial  (i.e.,  perpendi¬ 
cular  to  the  y-axis)  plane  in  the  uppermost  branches  of  the 
peripheral  airways.  The  velocity  fields  at  the  coronal  plane  show 
the  airflow  patterns  from  the  first  branch  (i.e.,  the  8th  generation)  to 
the  second  branch  (i.e.,  the  9th  generation)  in  the  peripheral  airways. 
In  the  normal  airway  model,  our  velocity  profile  was  skewed  toward 
the  inner  side  walls  in  the  second  branch  (Fig.  3(A)),  which  has  also 
been  observed  experimentally  in  previous  studies  [40,41],  In  the 
symmetric  obstruction  model,  the  flow  upstream  of  the  constricted 
zone  was  slowed  down  due  to  high  downstream  resistance  (Fig.  3 
(B)).  However,  while  passing  the  neck  of  the  constricted  zone,  the 
flow  evolved  into  a  jet-like  pattern  with  highly  increased  velocities, 
which  is  consistent  with  previous  studies  [14,21],  In  the  asymmetric 
obstruction  model,  the  airflow  was  disproportionate  with  ~  80%  of 
air  flowing  into  the  unobstructed  branch  (Fig.  3(C)).  In  the  random 
obstruction  model,  flows  in  the  second  branches  (i.e.,  the  9th 
generation)  in  the  peripheral  airways  were  unequal  even  though 
their  ainway  structures  were  the  same  (Fig.  3(D)).  The  flow  rate  in 
one  branch  was  50%  higher  than  the  other.  This  was  attributed  to  the 
uneven  distribution  of  obstructions  in  the  branches  downstream  of 
the  second  branch  (i.e.,  the  9th  generation)  (Fig.  2(E)).  With  the  same 
volume  reduction,  the  asymmetric  and  random  obstruction  models 
had  different  asymmetric  airflow  patterns.  This  suggests  that  the 
flow  pattern  is  not  only  dependent  on  the  degree  of  obstruction  but 
also  highly  dependent  on  their  distribution. 

We  compared  velocity  patterns  at  axial  planes  (i.e.,  along  the 
b-b')  in  the  normal  and  obstructed  airways.  Velocities  were 
highest  (i.e.,  colored  red)  for  the  normal  airway  model  and  lowest 
(i.e.,  colored  blue)  for  the  symmetric  obstruction  model.  The 
results  shown  in  Fig.  3  demonstrate  that,  although  the  velocity 
magnitudes  were  different,  the  velocity  patterns  for  steady-state 
inspiratory  flow  at  the  axial  planes  were  not  very  different 
between  the  models.  The  shapes  of  the  velocity  contours  in  all 
the  models  were  similar,  coaxial  circles,  with  the  highest  velocity 
at  the  center  and  zero  velocity  at  the  walls.  This  was  consistently 
observed  in  all  the  inspiratory  flows  driven  by  different  pressure 
drops  (data  not  shown). 

Fig.  4  shows  pressure  and  velocity  fields  for  high  expiratory 
flow,  driven  by  a  total  pressure  drop  of  18  Pa.  Similar  to  the  high 
inspiratory  flow  (Fig.  3),  a  higher  absolute  pressure  drop  was 
observed  for  the  obstructed  airway  models  than  the  normal 
airway  model,  reflecting  higher  resistance  in  the  obstructed  air¬ 
ways.  The  velocity  contours  in  the  coronal  plane  illustrated 
differences  in  the  upstream  flow  in  the  normal  and  obstructed 
airway  models.  In  the  normal  airway  model,  the  flow  merged 
symmetrically  and  gradually  from  the  lower  branches  after  pas¬ 
sing  the  branching  point.  This  resulted  in  a  peak  velocity  (colored 
red)  of  the  flow  being  increasingly  concentrated  around  the 
centerline  as  the  flow  progressed  toward  the  trachea.  This  has 
been  consistently  observed  in  experimental  studies  [40,42].  The 
increased  velocity  in  the  upper  branch  was  due  to  a  decreased 
total  cross-sectional  surface  area  [20,32].  The  flow  in  the 
obstructed  airway  models  was  not  as  smooth  and  gradual  as  in 
the  normal  airway  model.  Similar  to  the  inspiratory  flow.  Jets  were 
formed  in  the  neck  of  the  constricted  zones,  as  consistently 
predicted  in  other  computational  studies  [14,43],  Recirculating 
flows  were  present  in  the  vicinity  of  the  Jet  due  to  boundary  layer 
separation  (Fig.  4).  Vorticity  is  a  measure  of  the  rotation  of  the 
fluid  and  is  given  by  the  curl  of  the  velocity  field.  The  vorticity  in 
the  constricted  branches  of  the  symmetric  obstruction  model  was 
~  100-fold  higher  than  the  value  at  the  same  location  in  the 
normal  airway  model.  We  attributed  the  increase  in  vorticity  in 
the  symmetric  obstruction  model  to  the  significant  increase  in  the 
velocity  magnitude  near  the  constriction,  as  discussed  above 
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(Fig.  4)  as  well  as  to  the  increase  in  secondary  flow  velocity 
components  due  to  the  abrupt  change  in  the  shape  of  the  conduit 
near  the  constriction.  In  the  symmetric  obstruction  model,  the 
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Fig.  5.  Velocity  fields  of  low  expiratory  steady  flow  in  the  normal  airway  (A), 
symmetric  obstruction  (B).  asymmetric  obstruction  (C),  and  random  obstruction 
(D)  models.  The  applied  pressure  drop  across  the  extended  geometry  was  2  Pa.  (For 
interpretation  of  the  references  to  color  in  this  figure  legend,  the  reader  is  referred 
to  the  web  version  of  this  article.) 


velocity  dropped  immediately  after  the  two  jets  merged  at  the 
branching  point.  In  the  asymmetric  model,  the  flow  was  dom¬ 
inantly  from  the  unobstructed  branch,  similar  to  the  inspiratory 
flow  (Fig.  3(C)).  The  flow  rate  in  the  unobstructed  branch  was  80% 
higher  than  in  the  obstructed  branch.  For  the  random  obstruction 
model,  the  flow  rate  in  the  right  branch  (with  14  downstream 
obstructions)  was  50%  higher  than  the  left  branch  (with  18 
downstream  obstructions).  The  random  obstruction  model  also 
showed  the  recirculation  zones  in  the  second  branch  (i.e.,  9th 
generation)  of  the  peripheral  ainways. 

In  contrast  to  the  inhalation  flow,  we  obseiT/ed  widely  different 
axial  velocity  patterns  between  tbe  normal  and  the  obstructed 
airway  models  for  expiratory  flow.  As  shown  in  Fig.  4(A)  the  high- 
velocity  contours  (i.e.,  colored  red)  were  dumbbell  shaped  across 
the  y-axis.  Four  vortices  (marked  by  “-i-”)  were  formed  symme¬ 
trically,  as  has  been  found  experimentally  [40],  The  flow  in  the 
symmetric  obstruction  model  also  had  four  vortices  in  the  axial 
plane  of  the  uppermost  branches  but  with  much  lower  velocity. 
High-velocity  flow  was  more  concentrated  along  the  dorsal  and 
ventral  walls  (Fig.  4(B)).  In  the  asymmetric  obstruction  model,  the 
high-velocity  domain  was  highly  skewed  toward  the  obstructed 
branch  (Fig.  4(C)).  This  was  because  of  tbe  unbalanced  high- 
velocity  Jet  from  the  unobstructed  branch.  In  the  random  obstruc¬ 
tion  model,  the  velocity  contours  were  asymmetric  in  both  coronal 
and  axial  planes,  reflecting  the  uneven  distribution  of  obstruction 
sites  in  the  geometry  (Fig.  4(D)).  The  high-velocity  region  was 
skewed  toward  one  side  of  the  wall  with  the  formation  of  two 
large  vortices.  Also,  we  found  that  axial  airflow  patterns  were 
similarly  distinguishable  between  normal  and  obstructed  airways 
in  other  planes  as  well  (data  not  shown). 

We  examined  the  dependence  of  flow  intensity  on  airway 
patterns.  Fig.  5  shows  velocity  contours  of  low  expiratory  flow 
driven  by  a  2-Pa  total  pressure  drop.  Similar  to  the  high  expiratory 
flow  (Fig.  4),  Jet  flows  and  recirculation  zones  were  present  at  the 
obstructed  branches  but  with  smaller  degrees.  However,  the 
velocity  contours  in  the  axial  plane  were  not  distinguishable 
between  the  models  except  for  their  magnitudes.  Peak  velocity 
was  observed  at  the  center  of  the  axial  plane  with  little  deviation 
in  the  different  models.  Vortices  were  not  found  in  the  axial  plane 
due  to  reduced  secondary  flow  momentum. 


Table  2 

Similarity  measures  based  on  Pearson’s  correlation  coefficient  against  normal  velocity  contour  plots  at  the 
axial  plane  in  the  uppermost  branches  for  expiratory  flow. 


Symmetric 

Obstruction 

Asymmetric 

Obstruction 

Random 

Obstruction 

No  filter  (r) 

18 

0.94 

0.90 

0.83 

10 

0.98 

0.84 

0.95 

2 

0.99 

1.00 

1.00 

Low-velocity  filter  (rd 

18 

0.47 

0.41 

0.27 

10 

0,87 

0.21 

0.59 

2 

0.95 

0.81 

0.97 

Near-wall  domain  filter  (r2) 

18 

0.77 

0.65 

0.55 

10 

0.95 

0.50 

0.84 

2 

1.00 

0.96 

1.00 
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3.3.  Quantitative  characterization  of  airflow  patterns 

We  observed  marked  differences  and  distinct  characteristics 
in  flow  patterns  between  the  models.  Velocity  contour  plots 
(Figs.  3-5)  represent  not  only  the  spatial  distribution  of  velocity 
magnitudes  but  also  illustrate  the  airflow  patterns,  which  are 
given  by  the  distribution  and  shape  of  the  velocity  contours,  e.g., 
the  velocity  pattern  for  expiration  in  the  normal  airway  in  Fig.  4 
was  dumbbell  shaped.  Although  similarities  and  differences  in 
airflow  patterns  can  be  directly  visualized  from  such  plots,  a  more 
quantitative  characterization  of  flow  patterns  would  provide  for  a 
more  accurate  and  systematic  comparative  analysis. 

To  quantitatively  compare  flow  characteristics  between  the 
obstructed  airway  flow  patterns  and  the  normal  ainway  patterns, 
we  calculated  their  Pearson’s  correlation  coefficient.  The  correlation 
coefficient  measures  the  linear  dependence  between  two  variables, 
attaining  values  between  - 1  and  1,  with  - 1  indicating  perfect 
anticorrelation,  0  indicating  no  correlation,  and  1  indicating  perfect 
correlation  [34],  If  two  sets  of  velocity  contours  have  similar 
patterns  except  for  their  magnitudes,  the  correlation  coefficient 
would  be  close  to  1.  If  two  flow  patterns  share  no  similarity,  they 
would  be  uncorrelated  and  the  correlation  coefficient  would  be 
close  to  zero. 

Table  2  shows  correlation  coefficients  of  the  velocity  contour 
plots  at  the  axial  plane  in  the  uppermost  branches  for  expiratory 
flows.  We  obtained  these  values  by  comparing  the  velocity 
contours  from  the  obstructed  airways  with  those  from  the  normal 
airways.  For  all  flow  rates,  the  correlation  coefficient,  r,  was  close 
to  1  in  all  the  obstructed  airway  models  even  though  the  velocity 
contours  exhibited  significant  visual  differences  in  flow  patterns 
between  the  models  (Fig.  4).  This  was  mainly  a  consequence  of  the 
no-slip  boundary  condition  imposed  on  the  walls,  which,  in  effect, 
generated  a  dominating  low-velocity  field  near  the  walls  for  the 
normal  as  well  as  the  obstructed  airway  models  (Figs.  3-5). 

To  avoid  the  dominating  effect  of  this  near-wall  velocity  field, 
we  filtered  the  low-velocity  contours  and  recalculated  correlation 
coefficients  between  airway  models.  Fig.  6  shows  the  velocity 
contours  for  high  expiratory  flow  ( AP=  18  Pa)  after  the  use  of  two 
different  filters:  a  low- velocity  filter  (A)  and  a  near- wall  filter  (B). 
The  low-velocity  filter  eliminated  velocity  contours  with  velocity 
values  <  50%  of  the  maximum  velocity  to  limit  the  focus  to  the 
high-velocity  region.  Because  the  low-velocity  filter  produced 
different  areas  in  the  normal  and  obstructed  geometries  (Fig.  6 
(A)),  we  used  only  the  intersecting  areas  to  calculate  the  correla¬ 
tion  coefficient.  The  near-wall  filter  excluded  the  velocity  values  in 
the  region  close  to  the  airway  wall.  Because  the  boundary  layer 
thickness  [44]  of  the  flow  in  our  study  was  ~  14%  of  the  airway 
radius,  we  used  15%  of  the  airway  radius  as  the  distance  threshold 
for  the  near-wall  filter. 

Correlation  coefficient  values  obtained  after  applying  the  low- 
velocity  filter  (ti)  and  the  near-wall  filter  (r2)  were  dependent  on 
the  flow  conditions  as  well  as  the  airway  models  (Table  2).  For  the 
symmetric  obstruction  model,  correlation  coefficients  were  higher 
than  those  for  the  other  obstructed  models  for  all  flow  conditions 
(Table  2).  This  was  attributed  to  the  symmetry  of  airway  geome¬ 
tries  in  both  the  normal  and  symmetric  obstruction  models,  which 
led  to  similar  airflow  patterns.  For  high  expiratory  flow 
(AP=18Pa),  both  r-i  and  r2  were  <0.8  in  all  obstructed  airway 
models  (Table  2).  We  observed  minimum  ri  and  r2  values  for  the 
random  obstruction  model,  while  the  symmetric  obstruction 
model  resulted  in  the  maximum  values  (Table  2).  For  the  asym¬ 
metric  and  random  obstruction  models,  both  r-i  and  r2  were 
<  0.65,  demonstrating  highly  dissimilar  flow  patterns  compared 
with  the  normal  airway.  Importantly,  the  airway  resistance  values 
of  these  two  obstruction  models  were  merely  20%  higher  than  the 
normal  airway  model  (Table  1),  which  is  within  the  normal 
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Fig.  6.  Velocity  contours  in  the  region  of  interest  (ROl)  for  the  calculation  of 
correlation  coefficients  after  low-velocity  (A)  and  near-wall  (B)  filtering.  Velocity 
contours  in  the  uppermost  axial  planes  for  high  expiratory  flow  (AP=18  Pa)  were 
normalized  by  the  maximum.  The  ROl  was  redefined  as  velocity  fields  larger  than 
50%  of  the  maximum  for  low-velocity  filtering  and  as  15%  away  from  the  wall  for 
near-wall  filtering.  Correlation  coefficients  after  low-velocity  filtering  (r,)  and  near¬ 
wall  filtering  (r2)  were  obtained  between  the  velocity  fields  from  the  normal  airway 
model  and  each  obstruction  model.  (For  interpretation  of  the  references  to  color  in 
this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  article.) 


variability  range  of  airway  resistance  values  [45].  This  implies 
that,  for  high  expiratory  flow  and  after  applying  appropriate  filters 
to  avoid  inclusion  of  low  and  near-wall  velocities,  airflow  patterns 
are  a  much  more  sensitive  indicator  of  airway  conditions  com¬ 
pared  with  resistance  values  in  the  peripheral  airways. 

Conversely,  for  low  expiratory  flow  (AP=2Pa),  correlation 
coefficient  values  of  ti  and  r2  were  close  to  1,  indicating  that  the 
flow  patterns  were  similar  to  the  normal  aiiTAiay.  These  results 
were  consistent  even  if  different  thresholds  were  applied  for  the 
two  filters  (data  not  shown).  This  suggests  that  flow  patterns  with 
low  flow  rates  would  be  indistinguishable  between  the  normal 
and  obstructed  airway  condition,  which  was  also  visually  observed 
(Fig.  5). 

For  inspiratoiy  flow,  the  correlation  coefficients  at  the  upper¬ 
most  axial  plane  remained  1  for  all  flow  rates  even  after  we  used 
the  filters.  We  compared  flow  patterns  in  other  axial  planes  of  the 
uppermost  branches  and  tested  different  thresholds  for  both 
filters,  but  obtained  the  same  results.  This  corroborates  our  visual 
observation  that  the  airflow  patterns  at  the  axial  plane  of  the 
uppermost  branch  for  inspiratory  flow  are  similar  across  the 
normal  and  obstructed  airways. 
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Fig.  7.  Wall  shear  stress  distribution  in  airway  walls  of  the  normal  (A)  and  (B)  and  the  random  (C)  and  (D)  obstruction  models  for  high  inspiratory  and  expiratory  flows 
(AP=18  Pa).  The  locations  of  maximum  shear  stress  are  shown  in  the  insets.  (For  interpretation  of  the  references  to  color  in  this  figure  legend,  the  reader  is  referred  to  the 
web  version  of  this  article.) 


3.4.  Airway  constriction  and  wall  shear  stress 

Wall  shear  stress  is  the  tangential  force  per  unit  area  exerted  on 
the  airway  wall  due  to  airflow.  Fig.  7  shows  the  wall  shear  stress 
distribution  in  the  normal  and  obstructed  peripheral  airways  for 
inspiratory  and  expiratory  flows  with  the  pressure  drop  of  18  Pa. 
For  the  normal  airway  model,  wall  shear  stress  was  higher  in  the 
vicinity  of  the  bifurcation  compared  to  the  tubular  zone  in  each 
branch  (Fig.  7(A)).  For  inspiratory  flow,  the  maximum  wall  shear 
stress  was  observed  at  the  edge  of  the  inner  wall  of  the  first  branch 
(i.e.,  the  8th  generation)  where  the  flow  separates  into  two 
daughter  branches.  Because  wall  shear  stress  is  proportional  to 
the  velocity  gradient  perpendicular  to  the  flow  direction,  the 
location  of  the  maximum  shear  stress  was  consistent  with  the 
skewed  high-velocity  region  observed  near  the  inner  wall  for 
inspiratory  flow  (Fig.  3(A)).  Moreover,  the  velocity  gradient  in  the 
flow  was  largest  at  the  uppermost  branch  corresponding  to  the 
highest  flow  rate.  For  expiratory  flow  (Fig.  7(B)),  higher  values  of 
wall  shear  stress  compared  with  other  regions  in  the  airway  were 
observed  along  the  center  of  the  upper  branch  where  the  airflows 
from  the  daughter  branches  merge  and  result  in  a  high-velocity 
region  (Fig.  4(A)).  The  maximum  value  of  wall  shear  stress  for  the 
normal  airway  model  during  exhalation  was  near  the  bifurcation 
between  the  first  branch  (i.e.,  the  8th  generation)  and  the  second 
branches  (i.e.,  the  9th  generations)  of  the  peripheral  airways 
corresponding  to  the  highest  flow  rate  region  in  the  airway. 

In  the  obstructed  airway  models,  contrary  to  the  normal 
airway,  the  wall  shear  stress  was  higher  at  the  ain/vay  walls 
constituting  the  constrictions  for  each  branch  compared  with  the 
bifurcations  for  both  inhalation  and  exhalation  (Figs.  7(C)  and  (D)). 
At  the  locally  constricted  passages,  jet-like  flows  developed  with 


increased  stream  velocity  in  the  obstructed  airways  (Figs.  3  and  4). 
This,  in  turn,  induced  higher  velocity  gradients  along  the  radial 
direction  compared  to  other  locations,  increasing  the  wall  shear 
stress.  In  all  obstructed  airway  models,  the  maximum  value  of 
shear  stress  was  obseiA/ed  on  the  constricted  airway  wall  at  the 
uppermost  branch,  i.e.,  at  the  9th  generation  for  the  symmetric 
and  asymmetric  obstruction  models  (results  not  shown)  and  at  the 
10th  generations  for  the  random  obstruction  model  (Fig.  7(D)). 

Table  3  shows  the  maximum  wall  shear  stress  in  each  aiiway 
model  for  different  rates  of  inspiratory  and  expiratory  flows.  We  found 
that  the  wall  shear  stress  was  dependent  on  the  flow  rate  for  the 
same  airway  condition:  the  maximum  wall  shear  stress  increased,  on 
average,  by  ~  14-fold  for  a  ~  6-fold  increase  in  flow  rate  (see  Table  1 ). 
For  inspiratory  flow,  the  maximum  shear  stress  values  for  the 
symmetric  obshTiction  model  were  higher  than  those  of  the  normal 
airway  model  due  to  the  increase  in  flow  velocity  near  the  constric¬ 
tion,  as  discussed  above.  Similarly,  in  the  asymmetric  and  random 
obstruction  models,  the  maximum  shear  stresses  were  also  induced 
near  the  constrictions  in  the  obstructed  branches.  However,  due  to  the 
presence  of  unobstructed  branches  in  these  models,  the  airflow  rates 
through  the  obstructed  branches  were  lower  compared  with  those  in 
the  symmetric  obstmction  model.  Therefore,  the  maximum  shear 
stress  values  for  the  asymmetric  and  random  obstmction  models  were 
lower  compared  with  the  symmetric  obstmction  model  (Table  3). 

For  expiratory  flow,  the  maximum  wall  shear  stress  values 
were  higher  for  all  the  obstructed  models  compared  with  the 
normal  airway  model.  Comparing  across  the  obstructed  models, 
the  maximum  wall  shear  stress  value  was  the  highest  for  the 
symmetric  obstruction  model,  although  the  corresponding  flow 
rate  was  the  lowest.  This  was  again  due  to  the  disproportionate 
distribution  of  airflow  through  the  obstructed  and  unobstructed 
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Table  3 

Maximum  wall  shear  stress  in  the  normal  and  obstructed  airways. 


Maximum  wall  shear  stress  (Pa) 


Normal 

Symmetric 

Obstruction 

Asymmetric 

Obstruction 

Random 

Obstruction 

Inspiratory  flow, 

AP  (Pa) 

18 

0.49 

0.55 

0.23 

0.44 

10 

0.24 

0.32 

0.11 

0.22 

2 

0.03 

0.04 

0.02 

0.04 

Expiratory  flow, 

AP  (Pa) 

18 

0.21 

0.52 

0.26 

0.39 

10 

0.10 

0.27 

0.13 

0.21 

2 

0.01 

0.04 

0.02 

0.04 

blanches  in  the  random  and  asymmetric  obstruction  models. 
Between  the  random  and  the  asymmetric  obstruction  models 
whose  flow  rates  were  similar,  the  maximum  wall  shear  stress 
value  was  higher  for  the  random  obstruction  model. 


4.  Discussion 

4.1.  Airflow  patterns  and  airway  resistance  for  normal  and 
obstructed  airways 

In  this  study,  we  developed  normal  and  obstructed  CFD  models 
of  the  peripheral  lung  region  spanning  the  8th  to  the  14th  airway 
branches.  To  understand  the  pathophysiology  of  lower  ain/vay 
disease  conditions,  we  ascertained  how  the  airflow  patterns  were 
perturbed  as  a  function  of  airway  obstructions  under  different 
inspiratory  and  expiratory  flow  rates.  We  found  that  the  airflow 
velocity  and  pressure  magnitudes  were  significantly  different 
between  the  normal  and  obstructed  airways,  as  expected,  due 
to  the  increased  resistance  to  flow  in  the  obstructed  airways 
(Figs.  3  and  4).  More  importantly,  the  airflow  patterns  (the  shape 
of  velocity  contours)  were  found  to  be  markedly  different  both 
qualitatively  and  quantitatively  for  flow  rates  corresponding  to 
more  rigorous  breathing  conditions  compared  with  the  resting 
condition.  The  most  important  finding  from  our  work  was  that  the 
flow  profiles  demonstrated  clear  differences  between  normal  and 
various  disease  conditions  during  exhalation  (Fig.  4).  Moreover, 
obstructions  in  relatively  lower-generation  branches  were  mani¬ 
fested  in  airflow  patterns  in  the  corresponding  upper  ain/vay 
branches.  These  results  not  only  help  the  understanding  of  airflow 
characteristics  in  lower  airway  diseases  but  could  also  potentially 
aid  in  developing  and  improving  diagnostic  techniques  [46,47] 
based  on  imaging  of  airflow  in  the  lungs. 

In  our  comparison  of  airway  resistance  and  airflow  patterns 
between  normal  and  obstructed  airways,  we  showed  that  after  the 
application  of  appropriate  filters  to  avoid  the  inclusion  of  low-  and 
near-wall  velocities,  airflow  patterns  had  higher  sensitivity  than 
aiiTA/ay  resistance  in  differentiating  ainway  conditions.  Although 
resistance  values  for  the  symmetrically  obstructed  airways  were 
up  to  80%  higher  than  those  for  the  normal  airways,  for  the 
asymmetric  and  random  obstruction  models,  the  airway  resis¬ 
tances  were  only  20%  higher  (Table  1).  For  the  asymmetric  and 
random  models,  which  reflect  more  realistic  cases,  the  resistance 
values  themselves  were  within  the  normal  variation  for  different 
breathing  conditions  [45].  The  resistance  we  computed  here 


represents  airway  resistance  for  only  a  portion  of  the  peripheral 
lung  region,  which  contributes  <  20%  to  the  total  airway  resis¬ 
tance  for  the  entire  airways  [48],  Hence,  we  do  not  expect  that  the 
types  of  obstruction  cases  modeled  in  our  study  would  be 
detectable  by  measurements  of  total  airway  resistance.  Conversely, 
airflow  patterns  between  normal  and  obstructed  airways  were 
distinctively  different  even  for  similar  resistance  values.  The 
asymmetric  and  random  obstruction  models  showed  more  inho¬ 
mogeneous  airflow  patterns,  distributing  jet-like  flows,  and  recir¬ 
culation  zones  in  local  regions  in  the  lungs  (Figs.  3-5)  and  were 
distinct  from  normal  geometries  as  gauged  by  correlation  coeffi¬ 
cient  values  close  to  0.5  for  high  expiratory  flows  (Table  2). 

One  limitation  of  this  work  is  that  although  histological 
measurements  representative  of  obstructive  airway  conditions 
were  included  in  the  model  in  terms  of  a  reduced  surface  area 
[26,28],  we  could  not  include  the  architecture  of  the  obstructions 
due  to  unavailability  of  clinical  data.  We,  therefore,  used  idealized 
geometries  to  represent  the  obstructions  (Fig.  2).  However, 
because  we  used  smooth  transitions  of  the  geometry  from  the 
unobstructed  airway  sites  to  the  obstructions,  the  effect  of  the 
obstruction  architecture  on  the  flow  profiles  was  minimized.  The 
actual  obstructions  are  expected  to  be  irregular  and  would,  there¬ 
fore,  demonstrate  even  more  heterogeneous  changes  in  flow 
patterns  compared  to  the  normal  unobstructed  airway.  With 
advancement  in  imaging  methods,  it  may  be  possible  to  measure 
the  actual  structure  of  the  airway  obstructions,  which  could  then 
be  incorporated  in  the  model.  Another  limitation  of  the  work  was 
that  we  assumed  that  both  the  healthy  and  diseased  lungs  would 
experience  the  same  pressure  differentials  for  the  same  breathing 
effort,  which  may  not  be  valid  for  severely  diseased  cases  [49,50], 
Although  the  pressure  differential  needed  to  achieve  adequate 
alveolar  ventilation  should  be  higher  for  the  COPD  lungs  than  for 
healthy  cases,  the  capacity  of  the  respiratory  muscles  to  generate 
pressure  decreases  due  to  various  biomechanical  factors  in  dis¬ 
eased  subjects  [51  ].  As  a  result,  it  is  difficult  to  correlate  the  change 
in  pressure  differential  with  the  severity  of  the  disease.  In  the 
absence  of  such  data,  we  attempted  to  characterize  the  change  in 
airflow  characteristics  using  the  same  pressure  differential  for 
both  healthy  and  diseased  lungs. 

4.2.  Wall  shear  stress  generated  by  airflow  in  obstructed  airways 

The  wall  shear  stresses  generated  by  airflow  in  the  respiratory 
airways  induce  biologically  relevant  signals,  such  as  a  release  of 
adenosine  triphosphate  (ATP)  and  an  increase  in  intracellular 
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calcium  concentration  [52,53],  that  control  the  airway  defense 
mechanisms.  Moreover,  wall  shear  stresses  at  low  levels  may 
enhance  epithelial  barrier  function  [52],  In  contrast,  excessive  wall 
shear  stresses  may  result  in  the  damage  of  epithelial  cells  during 
airway  collapse  and  reopening  [54,55],  Therefore,  it  is  essential  to 
quantify  the  wall  shear  stresses  induced  by  airflow  for  various 
disease  conditions  to  assess  the  effect  of  obstruction  on  the  airway 
physiology.  In  the  absence  of  experimental  methods  to  estimate 
wall  shear  stresses  during  respiration,  computational  modeling 
provides  a  critical  recourse.  Although  modeling  efforts  have  been 
made  to  estimate  the  wall  shear  stress  in  the  central  airways 
[56,57],  it  has  not  been  characterized  for  iower  airways.  In  this 
study,  we  calcuiated  shear  stress  induced  in  the  airway  walls  in 
the  8th  to  the  14th  generations  as  a  function  of  breathing  rate  and 
distribution  of  airway  obstruction. 

We  observed  that,  for  normal  ainways,  the  wall  shear  stress  was 
highest  near  the  bifurcations  for  both  inspiratory  and  expiratory 
flows,  as  consistently  shown  in  previous  studies  [56,57].  However, 
for  obstructed  airways,  the  maximum  wall  shear  stress  was 
induced  near  the  constrictions.  Consequently,  one  would  expect 
the  increased  mechanical  forces  to  result  in  aggravated  patholo¬ 
gical  conditions  in  the  regions  of  constrictions  compared  to 
unobstructed  airways,  as  proposed  in  the  past  [58].  However,  our 
model  demonstrated  that  for  flow  rates  <  2300  ml/s  (i.e.,  flow 
rates  corresponding  to  heavy  activities),  the  wall  shear  stress  in 
both  normal  and  obstructed  airways  was  less  than  0.3  Pa,  which  is 
within  the  physiological  limit  needed  to  promote  respiratory  defense 
mechanisms  by  enhancing  epithelial  barrier  function  [52[.  Xia 
et  al.  [56]  showed  that  wall  shear  stress  in  compliant  airways 
decreases  by  50%  compared  with  rigid  ain/vays.  Therefore,  even  for 
flow  rates  >  4000  ml/s  (corresponding  to  forced  respiration  rates), 
the  maximum  wall  shear  stress  of  ~  0.5  Pa  calculated  by  our  model 
would  be  reduced  if  we  considered  airway  compliance,  which  was 
neglected  in  our  model.  This  implies  that  even  during  heavy  activities, 
airflow-induced  wall  shear  stress  would  not  elicit  damage  to  the 
airways,  but  rather  would  improve  defense  mechanisms  by  promot¬ 
ing  epithelial  barrier  function  [52[.  To  develop  a  better  understanding 
of  airway  injury  induced  by  shear  stress,  it  is  important  to  incorporate 
the  liquid-solid  interaction  in  the  airway  wall  surfaces  as  well  as  the 
airway  compliance.  It  has  been  suggested  that  epithelial  cell  damage 
occurs  due  to  increased  surface  tension  during  airway  reopening 
[54,55].  To  model  such  a  system,  it  would  be  necessary  to  develop  a 
multi-scale  model  to  describe  cellular  responses  in  addition  to  the 
global  airflow  distribution.  Our  work  provides  a  starting  framework 
to  develop  such  a  multi-scale  model. 

We  further  compared  the  airflow-induced  wall  shear  stresses 
for  different  distributions  of  airway  obstructions.  Although  the 
airway  diameter  was  reduced  by  the  same  amount  (50%)  at  the 
constricted  locations  for  all  obstructed  models,  the  maximum 
shear  stresses  were  widely  different.  In  the  symmetric  obstruction 
model,  the  wall  shear  stress  was  highest  compared  to  other 
obstructed  models  even  though  the  flow  rate  was  the  lowest. 
The  wall  shear  stress  in  the  random  obstruction  model  was  as 
much  as  2-fold  higher  than  the  asymmetric  model,  although  their 
flow  rates  and  airway  resistance  values  were  similar.  These  results 
demonstrate  that  the  distribution  of  obstruction  sites  is  an 
important  factor  in  determining  the  wall  shear  stress.  The  local 
anatomical  constrictions  and  the  interconnectivity  between  unob¬ 
structed  and  obstructed  branches  significantly  affect  the  accurate 
estimation  of  wall  shear  stresses,  making  the  use  of  3-D  models 
critical.  Although  earlier  models  could  describe  the  respiratory 
system  by  simplifying  the  3-D  lung  airways  into  a  network  of  one¬ 
dimensional  branches  [58,59],  the  details  of  3-D  flow  and  shear 
stress  distributions  as  presented  here  could  not  be  captured.  In 
this  study,  we  examined  3-D  characteristics  of  airflow  patterns  and 
shear  stress  on  the  airway  walls  by  comparing  normal  and 


obstructed  airway  models  from  the  8th  to  the  14th  generations. 
These  results  facilitate  the  understanding  of  how  lower  airway 
diseases  affect  inspiratory  and  expiratory  flow  phases  and  wall 
shear  stresses  for  various  distributions  of  airway  obstructions. 

5.  Summary 

In  this  study,  we  investigated  airflow  characteristics  in  the  lower 
always  of  healthy  and  diseased  human  lungs  using  computational 
models.  First,  we  developed  3-D  aiway  models  for  the  central  and 
peripheral  airways  with  different  numbers  of  airway  branches  and 
validated  the  computational  approach  of  using  reduced  number  of 
aiway  generations  to  model  respiratory  airflow  by  comparing  our 
results  with  in  vitro  experimental  data  from  the  literature.  Next,  we 
developed  3-D  computational  models  for  normal  and  obstmcted 
lower  airways.  We  included  symmetric,  asymmetric,  and  random 
obstructions  in  the  airways  between  the  8th  and  14th  generations 
and  compared  airflow  patterns,  airway  resistance,  and  wall  shear 
stress  under  a  range  of  breathing  conditions.  To  facilitate  a  quantita¬ 
tive  comparison  of  airflow  patterns  between  the  normal  and 
obstructed  airways,  we  obtained  the  Pearson's  correlation  coefficients 
for  cross-sectional  velocity  contours  after  applying  different  filters  to 
eliminate  the  low-velocity  field  near  the  airway  walls.  We  found  that, 
at  high  expiratory  flows,  the  obstmctions  in  relatively  lower- 
generation  branches  were  manifested  in  airflow  patterns  in  the 
corresponding  upper  airway  branches.  We  also  demonstrated  that 
airflow  patterns  subjected  to  filtering  exhibited  higher  sensitivity 
than  ain/vay  resistance  for  differentiating  obstructed  airways  from  the 
normal.  Further,  we  showed  that  wall  shear  stresses  were  not  only 
dependent  on  breathing  rates  and  degree  of  obstruction,  but  also  on 
the  3-D  distribution  of  obstmctions  in  the  lower  airways.  For  the 
same  degree  of  obstruction  and  breathing  rates,  we  noticed  as  much 
as  two-fold  differences  in  shear  stresses,  which  could  not  be  observed 
in  previous  one  dimensional  analytical  models.  In  contrast  to  previous 
studies  that  propose  increased  wall  shear  stress  due  to  obstmctions  as 
a  possibie  damage  mechanism  for  lower  airways,  our  models  demon¬ 
strated  that  for  flow  rates  corresponding  to  heavy  activities 
(  <  2500  ml/s),  the  wall  shear  stress  in  both  normal  and  obstructed 
always  was  less  than  0.3  Pa,  which  is  within  the  physiological  limit 
needed  to  promote  respiratory  defense  mechanisms.  In  conclusion, 
these  results  facilitate  an  understanding  of  how  lower  airway  diseases 
affect  inspiratory  and  expiratory  flow  phases  and  wall  shear  stresses. 
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